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F. Landis Markley* 

The quaternion has the lowest dimensionality possible tor a globally nonsingular attitude 
representation. The quaternion must obey a unit norm constraint, though, which has led to the 
development of an extended Kalman filter using a quaternion for the global attitude estimate and a 
three-component representation for attitude errors. We consider various attitude error representations 
for this Multiplicative Extended Kalman Filter and its second-order extension. 

INTRODUCTION 

Use of the Extended Kalman Filter (EKF) 1 ^ for spacecraft attitude estimation has a long history. The 
earliest application*^ used three Euler angles to parameterize the attitude, and one of the most recent 
examples employs the rotation vector. 6 Both of these parameterizations require transcendental function 
evaluations and are singular for certain attitudes, so other parameterizations have been employed for Kalman 
filtering. The search for an ideal parameterization is frustrated by the fact that a globally nonsingular three- 
dimensional parameterization of the rotation group is topologically impossible. 7 

The unit quaternion has the lowest dimension of any globally nonsingular attitude parameterization and 
represents the attitude matrix as a homogeneous quadratic function. The kinematic propagation equation for 
the quaternion is bilinear in the components of the quaternion and the angular velocity vector. These 
advantages have motivated the development of filters using the quaternion representation for both the global 
attitude and the attitude errors. 8 The linear measurement updates provided by these straightforward EKFs 
violate the nonlinear quaternion normalization constraint, however. Several solutions have been proposed 
for this problem, with varying degrees of success. 8-1 ^ 

In this paper, we will consider an EKF that parameterizes the global attitude with a unit quaternion, while 
employing a three-component representation for the attitude errors. 9,114 Several justifications have been 
provided for this non-standard approach to Kalman filtering. The approach in Section IX of Ref. 9 regards 
the filter as estimating a four-component quaternion, but argues that the quaternion covariance is rank- 
deficient due to quaternion norm constraint, or equivalently due to the underlying three-dimensional nature 
of the rotation group. 71116 Since the 4x4 quaternion covariance is assumed to have rank three, it can be 
projected onto a 3x3 matrix without any loss of information. This gives a measurement update obeying the 
quaternion norm constraint to first order in the measurement residual, so enforcing exact quaternion 
normalization only modifies the quaternion to second order and is therefore outside the purview of the EKF. 
Vathsal 17 has extended this approach to second order; but the result does not appear to be a consistent 
second-order filter, since his quaternion update violates the norm constraint in second order. 

An alternative approach, which has become known as the multiplicative EKF (MEKF), is set forth in 
Section XI of Ref. 9. This regards the filter as estimating a three-component attitude error vector, which 
naturally has a 3x3 covariance, while parameterizing the global attitude by a quaternion. This is actually the 
original approach, 1114 and it leads to the same algorithm as the covariance projection approach in the 
linearized EKF approximation. It rests on a firmer conceptual basis, however. 

This paper is an extension of Ref. 9, which provides citations of the literature through 1981. Its principal 
aims are to consider in some depth the relationship between the three-component attitude error vector and 
the quaternion, and to extend the MEKF approach to a consistent second-order Kalman filter. To set the 
stage, we begin with a brief discussion of attitude representations. We then develop the basic equations of 
the MEKF, and consider its extension to second order. 
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ATTITUDE PARAMETERIZATIONS 


Good reviews of attitude representations are available. 1819 We provide a brief discussion to establish 
conventions and notation for the representations of interest to this paper. We regard the 3x3 orthogonal 
attitude matrix, or direction cosine matrix as the most fundamental representation of the attitude. 


Rotation Vector 

Euler’s Theorem 20 states that the most general motion of a rigid body with one point fixed is a rotation by 
an angle 0 about some axis. We specify the axis by a unit vector e and combine the axis and angle into a 
rotation vector ♦ Me. The attitude matrix as a function of § is 

A(<|>) = / 3x3 -(sin0/0)[<|>x] + [(1 -cos0)/0 2 ][<|>x] 2 , (1) 


where the 3x3 identity matrix is denoted by / 3x3 and the cross product matrix is 


qx] 


0 

~<h 

<72 


0 

"<7. 


<7i 

0 


( 2 ) 


Equation (1) and the kinematic equation for the rotation vector are transcendental and ill behaved, though 
finite, for zero rotation angle. All rotations can be mapped to points inside and on the surface of a sphere of 
radius 7 iin rotation vector space, where points at opposite ends of a diameter represent the same rotation; 
but the rotation vector may jump from one end of a diameter to the other as the attitude varies smoothly. 
These jumps can be postponed by expanding the representation to a sphere of radius 2k, but they cannot be 
avoided entirely, since the kinematic equation for the rotation vector is singular for 0 = 2 it. These 
characteristics limit the usefulness of the rotation vector as a global attitude representation. 


Quaternions 


A unit quaternion representing spacecraft attitude has a three- vector part and a scalar part, which are related 
to the axis and angle of rotation by 



Qv 


q 


esin(0/2) 

<7 = 

_ c u_ 


_<?4 _ 


cos(0 / 2) 


and obey the unit length constraint 

M 2 = |qf + ?4 = i- 


(4) 


The four components of q are the Euler symmetric parameters or the Euler-Rodrigues parameters, which 
first appeared in a paper by Euler 21 and in unpublished notes by Gauss. 22 Rodrigues’ classic paper of 1840 
first demonstrated their general usefulness. 23 Hamilton introduced the quaternion as an abstract mathematical 
object in 1844, 24 but there is some question as to whether he correctly understood its relation to rotations. 25 

Unit quaternions reside on the three-dimensional sphere S3 embedded in four-dimensional Euclidean space 
E4. The attitude matrix a homogeneous quadratic function of the components of a unit quaternion; 

A (q) = (q\ - |qf )/ 3x3 - 2<7 4 [q x ] + 2qq r . (5) 


The quaternion representation is 2-1 because Eq. (5) shows that q and —q represent the same rotation matrix. 
The quaternion is an ideal global attitude representation, since it varies continuously over S3 as the attitude 
changes, avoiding the jumps required by some three-dimensional parameterizations. It is customary to 
restrict a quaternion representing an attitude error to the hemisphere of S3 with q A > 0, however. 
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We follow Ref. 19 in writing the quaternion product as 

^ >4q + <?4P-P X <ll 

p® q = ■ W 

L p 4 ^ - p q J 

This differs from the historical multiplication convention, 18 ' 24 denoted by qp without an infix operator, by 
the sign of the cross product in the vector part. The two products are related by p®q = qp . The convention 
adopted here has the useful property that 

A(p)A(q) = A(p®q), ( 7 ) 

as contrasted with A(p)A(q) = A(qp ) . Equation (7) means that the rotation group and the quaternion group 
are almost isomorphic, the qualifier “almost” owing to the 2-1 nature of the mapping. 16 We use an overbar 
to denote the quaternion representation of a three-vector. 

(8) 


With this convention, the kinematic equation for the quaternion can be written in the alternative forms 

q = 7(0® q = ^£2((0)<?, (9) 

where G) is the angular velocity vector, and the skew-symmetric 4x4 matrix Q(to) is defined by 


Q(a>) s 


-[<ox] 

-co r 


(0 

0 


(10) 


The skew-symmetry of Q(a>) preserves the normalization of q , but this normalization may be lost due to 
computational errors. If so, it can be restored trivially by q = q/\q\. 


Gibbs Vector or Rodrigues Parameters 


The three Rodrigues parameters are defined by 23 


q _ esin(0/2) 
q A cos(0/2) 


e tan(0/2) . 


(ID 


Gibbs arrayed them in a vector semitangent of version', it is little wonder that we now call it the Gibbs 
vector. 26 The Gibbs vector can be regarded as a gnomonic projecton of the S3 quaternion space onto three- 
dimensional Euclidean g space, as shown in Fig. 1. This is a 2-1 mapping of S3, with q and -q mapping 
to the same point. Since q and — q represent the same rotation, the Gibbs vector parameterization is a 1—1 
representation of the rotations onto E3. The Gibbs vector is infinite for 180° rotations (the q 4 = 0 equator 

of S3), which is undesirable for a global representation of rotations. The attitude matrix has the Gibbs 
vector representation 

^(g) = /3x3- 2 ( 1 + « 2 )"'{tgX]-[g><] 2 }. ( 12 > 


where italics are used to denote the magnitude of any three-vector other than the vector part of a quaternion. 


Modified Rodrigues Parameters 


These parameters were introduced by Wiener: 27 


P = 


q 


esin(0/2) 

1 + cos(0/2) 


etan(0/4) , 


(13) 
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(14) 


and represent the attitude matrix by 

A(p) = / 3x3 -4(1 + /? 2 ) _2 {(l-p 2 )[px] - 2 [px] 2 } . 

Marandi and Modi 28 pointed out that the Modified Rodrigues Parameters can be viewed as a stereographic 
projection of S3 quaternion space, as shown in Fig. 2. One hemisphere of S3 projects to the interior of the 
unit sphere in three-dimensional p space, and the other hemisphere of S3 projects to the exterior of the unit 
p-sphere. All rotations can be represented by Modified Rodrigues Parameters inside and on the surface of the 
unit ball. If we extend the representation to all Euclidean p space, we have a 2-1 parameterization with p 
and -p/p 1 representing the same rotation. This parameterization shares many characteristics with the 
rotation vector parameterization, including the need for discrete jumps, but avoids transcendental functions. 


g 


▲ </i 


Figure L Gibbs Vector 
as a Gnomonic Projection 



Figure 2. Modified Rodrigues Parameters 
as a Stereographic Projection 



MULTIPLICATIVE EXTENDED KALMAN FILTER 

The MEKF represents ihe true nonlinear state as the quaternion product 

q(t) = Sq(a(t))®q rt , f (t), (15) 

where q re/ (t) is some unit reference quaternion and &7(a(0) is a unit quaternion representing the rotation 
from q ref (t) to the true attitude parameterized by q{t). We parameterize 5q(a(t)) by a three-vector a (t) of 
attitude errors in the spacecraft body frame. The two attitude representations a(r) and q ref {t ) in Eq. (15) are 
clearly redundant. The MEKF computes an unconstrained estimate of the three-component a(f) while using 
the four-component q r f (t) to provide a globally nonsingular attitude representation. An alternative 
formulation, which has some advantages, reverses the order of multiplication in Eq. (15) so that a(f) 
represents the attitude errors in the inertial reference frame rather than the errors in the body frame. 29 

Given an estimate of the error vector a(r), the best estimate of the attitude quaternion is 

q(t) = Sq{k{t)) ® q r ef( f )' < 16 ) 

which is a unit quaternion. We remove the redundancy in the attitude representation by choosing the 
reference quaternion q rt>J ( t ) so that the expectation value a(r) is identically zero when all available dynamic 


4 



and measurement information is included. Then Eq. (16) shows that q re j (0 is identically equal to r/(r), so 
we will denote it as such. This means in turn that Sq(a(t)) is a representation of the attitude error. It is 
important to note that if q(t) had been defined as the expectation of some probability distribution over S3, 
it would not be a unit quaternion unless the distribution were concentrated at a point, since restricting the 
probability distribution in quaternion space to the surface of S3 means that its expectation value must be 
inside S3. This is the conceptual advantage of the MEKF. 

Continuous/discrete filtering proceeds in three steps: time propagation, measurement update, and reset. The 
continuous time propagation is arranged to keep a (t) = 0, but the discrete measurement update assigns a 
finite post-update value a(+) to a. In order to avoid the need to propagate two representations of the 
attitude, the reset operation moves the attitude information from a(+) to q(+), after which a is reset to 
zero. Since true quaternion is not changed by this operation, Eq. (15) requires 

&/(a(+))® q{-) = & 7 ( 0 )® <?(+) = <?(+), (17) 

where q(-) is the pre-update quaternion. It is possible eliminate the discrete reset by keeping a (f) = 0 at all 
times, even during the update, by considering each attitude measurement update to be spread out over an 
infinitesimal time interval, rather than being instantaneous. 30 This paper treats measurement updates as 
discrete rather than continuous, however. 

The significance of the reset is obscured in the standard EKF, which represents the true state X as the sum 
of the reference value X and a small error x, 

X = X + x . (18) 


The measurement processing produces an updated value of the error vector 

x (+) = x(-) + Ax. (19) 

The reset operation moves the update information from the error state to the estimate of the full state by 

X(+) = X(-) + x(+) -x(-) = X(-) + Ax, (20) 

which gives the appearance that the update is applied directly to the full state estimate. In fact, an implicit 
reset of the gyro biases is preformed in exactly this manner in the filter considered later in this paper. The 
reset of the attitude must be treated explicitly in the MEKF, however. 


Attitude Error Representations 


One possible parameterization of a is the rotation vector 0 e, which we denote a^ , so from Eq. (3), 


Sq( a,) 


(a,/a,)sin(a, 12) 
cos (a+ 12) 


(21a) 


This parameterization has the advantage that the covariance includes the angular variances in radians 2 , but it 
is numerically inconvenient owing to the appearance of trigonometric functions and a 0 /a (> , which is 
indeterminate at a 0 =0. 

We can retain this inteipretation of the covariance matrix by requiring a to be equal to the rotation vector in 
the small angle approximation. This leads to the second parameterization of a as twice the vector part of 
Sq, which is the parameterization used in Section XI of Ref. 9, except for the factor of two; 


&?(*<)= ^ 



(21b) 
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A third alternative is four times the vector of Modified Rodrigues Parameters, so from Eq. (13), 


&7(a„) = 


16 + a; 


8a r 

16 -a^ 


(21c) 


This parameterization has the computational advantage of not requiring any transcendental functions. 
A fourth parameterization is twice the Gibbs vector; from the inverse of Eq. (11), 

&/(a f ) = 


1 

V 

^4 + a 2 s 

2 _ 


(2 Id) 


The Gibbs vector form has the advantage that the reset can first define the unnormalized quaternion 


P = 





and then the updated unit quaternion is given by 

q(+) = p/\p\. 


( 22 ) 


(23) 


This form avoids any accumulation of numerical errors in the quaternion norm. 


These four definitions of a provide the same second-order approximations to the quaternion, 


and to the attitude matrix 


Sq( a) = 


a/2 

1 -a 2 /8_ ’ 


A(8q) — / 3 x 3 [ax] + 2 ( aa a ^ 3 x 3 )- 


(24) 

(25) 


Thus they are equivalent for an EKF, which uses a linear approximation; but they differ in third and higher 
orders in a. It is worthwhile to note that Eqs. (24) and (25) only hold to first order if the components of a 
are taken to be Euler angle rotations about three orthogonal axes, as in Refs. 5 and 13. Thus, an Euler angle 
parameterization will lead to the same EKF, but will give different results for a second order filter. 

The error representations are quite different for large errors, which may arise before the Kalman filter has 
converged. Consider the largest possible attitude error, which is a 180° rotation. This corresponds in the 
various three-component representations to a q - 2, - n, a p = 4 or infinite a r A filter update a(+) larger 

than this value does not give a sensible quaternion reset. The infinite limiting value for the Gibbs vector 
representation is another reason to prefer this representation. 

State and Covariance Propagation 


The filter dynamics are found by differentiating 


using Eq. (9) and the identity 


Sq(t) = q(t) ® q~' (t) , 
dq~' /dt = —q~' ®q®q~'. 


(26) 

(27) 


This follows, like the corresponding expression for the derivative of the matrix inverse, from the vanishing 
of the time derivative of q® q~' , which is the constant identity quaternion. This gives 

d(<5q(a))/dt = jio ® 5q(a) - Sq(a) ® q ® q~' . (28) 
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Taking the expectation value of this equation in the usual linear EKF approximation of ignoring 
correlations gives _ 

d(<5q(a))/df = +cb ® c5^(a) -Sq(k)® q® q ~ x . (29) 

Since &/(a) is the identity quaternion, which is a constant, we have 

0 = (30) 


or 

q(t) = \<a(t)®q(t) = i£2(<B(0)gW, ( 31 ) 

This quaternion propagation equation is the same as the equation derived by more conventional methods, 
but we have derived it from the requirement that a(f) be identically zero. 

Now we specialize to the case where a set of gyros is used to obtain angular rate information in place of 
models of the spacecraft dynamics. The Kalman filter estimates the spacecraft attitude and a three-vector of 
gyro drift biases. 9 ' 14 We employ Farrenkopf s gyro dynamics error model, 31 which means that we ignore the 
output noise for rate-inlegrating gyros. 32 This is an excellent approximation for navigation-grade gyros. 

The angular rate vector is given in terms of the gyro output vector d>(r) and gyro drift vector b(/) in 
spacecraft body coordinates by 

to(/) =m(t)-b(t)-r\ l (t), (32) 

where "p, (r) is a zero-mean white noise process. The gyro drift vector obeys 

b(f)=r] 2 (r), (33) 

where i) 2 (r) is an independent zero-mean white noise process. This clearly implies that 

b(?) = 0 (34) 


We consider the Gibbs vector parameterization for specificity. Solving Eq. (2 Id) for a, gives 

a g (t) = 2(8q) v /(8q) i , (35) 

where the subscript V and 4 denote the vector and scalar parts of the quaternion, as in Eq. (3). Then, with 
Eqs. (15) and (28), 


f(x(i),/)sa ( (i) = ]^o)(/)® 


4 H(o(r)< 




M') 

2 


® q(t) ® q (r) 




a,(0 


(36) 


®q{t)® q “'(/)[ a At), 


where the time dependence of the reference quaternion is implicitly included in the time argument of 
f(x(/),f). Inserting Eq. (31) gives, in the EKF approximation 

f(x(/),0 = -cb x a ? (r) + Am(r) Aa»(r) x a^ (/) + Aco(r) • a^r)] a g (r) , (37) 

where 

Ao)(r) = o)(r) — co(r) =o>(r)-[ci>(0-b(/)] = b(r) - b(/) - tj, (/) . (38) 
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We want to construct an EKF for the six-component state vector 


x(/) = 


a(f) 

b(r) ' 


(39) 


where we suppress the subscript on a in equations that hold for any three-dimensional parameterization. 
This state vector obeys the propagation equation 


a 1 


‘f(x(0,0 

b 


th(') 


(40) 


The expectation value of the state vector was shown above to be constant in the absence of measurements 
The covariance matrix can be partitioned into 3x3 submatrices as 

>„(/) P c {t) 


P(t) = E{(x-x)(x-x) r ) = 


PjU) P h (0 


We let di/da denote the matrix with elements 

[df/daljsdfi/dcij 

and define the matrices 

~E{df/3a} E{r)f/3b} 


F(t)> 


G(t): 


^3x3 ^3x3 


Ejaf/ari,} E{3f/3ri 2 } 

0 3x3 / 3x3 


0 ^3x3 


Then the time propagation of the covariance is given by 

P{t) = F(t)P(t)+ P{t)F T (t) + G{t)Q(t)G T (t)- P(t)H T (t)R-'(t)H(t)P(t), 

where 

fe,(o o 3x3 


< 2(0 = 


®3x3 Ql(0 


(41) 


(42) 


-[co(/)x] 

f 3x3 

(43) 

^3x3 

^3x3 

— f 3x3 

^3x3 

(44) 


(45) 


(46) 


with 

E{T]/(0Tl^(r')} -02,(0 for ij= 1,2. (47) 

The covariance propagation during the intervals between the attitude measurement updates is identical to 
Refs. 13 and 14, except for factors of | in some of the formulations. The Appendix shows that these 
equations also hold if Eq. (21b) is used in place of Eq. (21d). By implication, these relations are independent 
of the representation chosen for a, as long as Eq. (24) is satisfied. 

Measurement Model and Update 

A measurement is modeled as a an ^/-dimensional function h of a vector \ B in the spacecraft body frame, 
with w'hite noise added; 

z = h(v fl ) + white noise. (48) 
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The representation of \ B in the body frame is the mapping of its representation v, in the inertial reference 
frame by the attitude matrix: 

v B =A(q)v, =(/ 3 * 3 -[ax])v s , (49) 


where we have used Eqs. (7), (15), and the first-order approximation to Eq. (24), and where \ B = A(q)\ , is 
the vector in the body frame that is predicted by the reference quaternion. This gives, to first order in a. 



h(v s ) = h(v„) -0h/3v)|.Jax]v e = h+ H a a, 

(50) 

where 


h = h(v B ) 

(51) 

and 


H u = Oh/3v)| [v,x], 

' Y a 

(52) 

Thus the 

m x 6 measurement sensitivity matrix is 



H = [3h/3a dh/3b] = [//„ 0 mx3 ] , 

(53) 


since the attitude measurements are assumed not to depend explicitly on the gyro drifts. The Kalman gain 
matrix is given by 

K = [P a H P, (-)f H T a [ H a P a {-)H T a + A]"' , (54) 

where R is the covariance of the measurement white noise. The state update is given by 

x(+) = x(-) + K( z-z) = x(— ) + K\ z-h - // fl a(-)], (55) 


where z denotes the measured value and z is the value predicted from the pre-update state estimate. The 
covariance update is 

P(+) = P(-) -[/»„(-) P c H] T H T a [H a P a (-)Hl +R]~' H u [P a (~) />(-)]. (56) 

Reset 

The quaternion reset uses Eq. (17) with any of Eqs. (21); we have seen that the update is independent of the 
representation chosen. Resets can be performed after each measurement update, in which case the term 
H a a(-) in Eq. (55) is identically zero; but the reset is usually delayed until all the updates for a set of 
simultaneous measurements have been performed, for computational efficiency. It is imperative to perform 
a reset before beginning the next time propagation, however, to assure that x(r) is zero at the beginning of 
the propagation, and thus to avoid the necessity of propagating \(t) between measurements. The reset does 
not modify the covariance, since it neither increases nor decreases the total information content of the 
estimate; it merely moves this information from one part of the attitude representation to another. 

Quaternion Measurements 

Many modern star trackers track between 5 and to 50 stars simultaneously, match them to stars in an 
internal star catalogue, and compute their attitude as a inertially-referenced quaternion. 33 34 The 
computation also produces an estimate of the 3x3 covariance of the attitude error angles. 35 36 It is a simple 
matter to transform these quantities from the star tracker reference frame to the spacecraft frame to produce a 
quaternion “measurement” q and a 3x3 measurement covariance matrix R. The most convenient way to 
present the information in the measured quaternion to the Kalman filter is in terms of one of the three- 
dimensional parametenzations of the estimated attitude error 

Sq(a -z) = q® q~ x (-). (57) 
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The measurement model is simply 


h(a) = a , (58) 

so H a is the 3x3 identity matrix and R is the covariance of this error angle. Since h is identically zero from 
Eqs. (57) and (58), the state update is given by 

a(+) - a(-) + K a [z- a(-)] = (/ 3x3 - K a ) a(-) + K a ~z (59a) 

and 

b(+) = b(-) + K h [z- a(— )], (59b) 

where 

K a =PA-)[P a H + Rr l , (60a) 

and 

K b =P^-)[P a (-) + RY l . (60b) 

It is important to use the same three-dimensional parameterization in the measurement model, Eq. (57), as 
is used in the reset. For example, if 

3j (/) = 2[ ~q <S> q' (-)]„/[ q ® q' (-)] 4 , (61) 


as in Eq. (35), then Eq. (2 Id) should be used for the reset. Alternatively, if the measurement is modeled as 
twice the vector part of the error quaternion, the reset should use Eq. (21b) in place of Eq. (2 1 d). With this 
proviso, we see that when R « P a (-), so that K a goes to the 3x3 identity matrix, we have q{+) = q. 

SECOND ORDER FILTER 

Second-order terms in the nonlinear Kalman filter can become important when nonlinearities are significant 
relative to the measurement and process noise terms. According to Maybeck, 4 a first-order filter with bias 
correction terms obtains the essential benefit of a second order filter without the computational penalty of 
additional second moment calculations. This filter adds second-order corrections to the state propagation and 
measurement residual equations, but uses the EKF expressions for the covariance and gains. 


State and Covariance Propagation 


The state estimate is propagated by 


x(t) = 


f(x(r),r) + b p (r) 

0 


(62) 


where the propagation bias correction is given by 




11 j dX'dXj " 


(63) 


For the Gibbs vector parameterization, Eqs. (36), (41) and (63) give 

b p (0 = i P„ ( ® (0 — 2 [ qit) ® q '' (t )] v } + CO r (0 , 


where (o r has components 


(64) 

(65) 
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with {i,j, /t} a cyclic permutation of ( 1, 2, 3}. Including this second order term in Eq. (36) gives 


a/?) = H©(?)® 




a,(r) 


a s (f) 




® q(t) ® q-' (?) | + \ P a {©(0 - 2[ q(t ) ® q (?)] v } 




( 66 ) 


q(t)®q~'(t)\ a^?) + ©,.(?). 


The condition that a (/) and a (?) are equal to zero is 


[ «j(?) ®q-' (?)], = | { W (?) + [ / 3x3 + 1 P 0 (?)]"' 0), (0 } , 


(67) 


It is shown in the Appendix that the factor of [/ 3x3 + ■J’ 7^,(?)] 1 depends on the specific choice of the three- 
dimensional parameterization of the rotation. Since P a {t) and co r are both second order in the estimation 
errors, it is consistent with a second-order filter to replace [/ 3x3 + } P a 0)] by the identity matrix, giving 


q(t ) = ^ ( <0 (?) + CO , (?)] ® q(t) = ^ Q(© (?) + 0) r (?)) q(t) . 


(68) 


Time propagation between measurements is changed from the EKF by the addition of this second-order 
correction to the angular rate vector arising from the skew part of the covariance between the attitude errors 
and gyro drift bias errors. This is equivalent to the result obtained by Vathsal. 17 


Measurement Update and Reset 

In the second-order filter, the predicted measurement z has the form 

z = h + H a a(~) + b m (r), (69) 


with the measurement bias term given by 4 




3 2 h(x,r) 


a,*, F < n 


x=x(0 


Using Eq. (25) to expand Eq. (48) to second order in a gives 


h(v g ) = h(v a ) + 


3h(v) 


r)v 


{[v 6 x]a + ^(aa r ~ a 2 1 b} + 2 ^ 


a 2 h(v) 


fn dv > dv j 


(70) 


([VeXlaMtVgXla),.. (71) 


Inserting this into Eq. (70), using the symmetry of P a and the mixed second-order partial derivatives and the 
fact that the measurement does not depend explicitly on the gyro drift bias, gives 

b„,(?) =±0h/3v)|. [P a - (tr P a )/ w ]v fl + + X (3 2 h/av,dv y ) |_ ([v B x] T P a [v „x]) „ , (72) 

/J. I 


where tr denotes the matrix trace. This result differs from the measurement bias found by Vathsal, whose 
computation ignored the quaternion norm constraint. 

Now consider the special case that the attitude covariance P a is a multiple of the identity, P a = pj 3x3 for 
some scalar p a . This case is of interest since a Kalman filter is often initialized with a large covariance of 
this form. We want to be sure that an unrealistically large covariance won’t corrupt the update. In this case 


b ro (f> ^ /V -0V 3v )|v/* 

i'vH V " 


(73) 
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Now consider two different measurement models for a star tracker. The first is the usual model, 


where 


h(v) = 


“l/«3 

Lm 2 /«,J’ 


U = B v. 


(74) 

(75) 


is the unit vector to the star in the star tracker reference frame, which is rotated from the spacecraft body 
frame by the orthogonal transformation matrix B. The measurement sensitivity matrix is 


H 


a 


1 w 3 



0 


~ U 2 

-w, 


[ux]fl. 


(76) 


The first term in Eq. (73) for the second order measurement bias vanishes, and the second term gives 


The second measurement model is linear: 


b .n(0 = P fl (l + A 2 )h. 


h(v) = 


The measurement sensitivity matrix for this measurement model is 

"l 0 O' 

0 1 0 

In this case, the second term in Eq. (73) for b m (f) vanishes, and the first term gives 


H . = 


[uxjfi. 


(77) 

(78) 

(79) 


b m (0 = -p> (80) 

These two measurement models give measurement biases of the same order of magnitude but with opposite 
signs. For larger than usual star tracker initialization errors of 0.01 radians, or 0.573°, the correction to the 
predicted measurement is 0. 1 % of the leading term h. 

Since the input to the Kalman Filter for a quaternion measurement is linear in the three-dimensional attitude 
parameter vector, the measurement bias b m (t) is identically zero in this case. Also, since Eqs. (21) are 
exact to all orders of a, the reset for the second-order filter is the same as the EKF. 


SUMMARY 

The major result of this paper is to clarify the relationship between the four-component quaternion 
representation of attitude and the three-component representation of attitude errors in the widely used 
estimator that has become known as the Multiplicative Extended Kalman Filter. We view this filter as 
based on an apparently redundant representation of the attitude in terms of a reference quaternion and a three- 
vector specifying the deviation of the attitude from the reference. This apparent redundancy is removed by 
constraining the reference quaternion so that the expectation value of the three- vector of attitude deviations 
is identically zero. It is therefore not necessary to compute this identically zero expectation value. The basic 
structure of the Multiplicative Extended Kalman Filter follows from constraining the reference quaternion in 
this fashion: the reference quaternion becomes the attitude estimate, the three- vector becomes the attitude 
error vector, and the covariance of the three-vector becomes the attitude covariance. All these results are well 
known in practice, but the justification for using this mixed attitude representation has been unclear. 

Several different three-dimensional parameterizations give identical results in the linear Extended Kalman 
Filter and in a second-order filter, except in the reset step where they differ in third order in the measurement 
update. The second-order propagation and measurement biases are easily computed in this framework. 
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APPENDIX 


Instead of using the Gibbs vector, we will employ the vector part of the quaternion as in Eq. (21b). With 
this parameterization, Eq. (36) is replaced by 


f(x(/),r) = a,(/) = Hffl(t)® 


a „(0 


a, ; (f) 

ft -a)(i) 


®q(t)®q~' (r)| 


(Al) 


In the EKF approximation of Eq. (31), 

f(x(/),0 = -to x a q (t) + ’2-^4 - a* A©(/)-^Aa)(r)xa <p (r), (A2) 


which is the same as Eq. (37) through terms of second order in a and A© . Thus the EKF propagation and 
update equations are unchanged. Differentiating Eq. (Al) gives, with Eq. (63), 

b p (/) = - i trP B (/) ( © (, t ) - 2[ q(t) ®q~' (t )] v } + © r (/) , (A3) 

in place of Eq. (64), which changes Eq. (67) to 

[<7(0® <7 (f)] v = 2 {©(0 + [l - i tr P a (0]"'®, (0) • (A4) 

Equations (67) and (A4) agree if and only if we ignore products of P a (t) and © r . 
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